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ABSTRACT 

We present a maximum-entropy method (MEM) and 'Mexican Hat' wavelet 
(MHW) joint analysis to recover the different components of the microwave sky 
from simulated observations by the ESA Planck satellite in a small patch of the sky 
(12.8 x 12.8 deg 2 ). This combined method allows one to improve the CMB, Sunyaev- 
Zel'dovich and Galactic foregrounds separation achieved by the MEM technique alone. 
In particular, the reconstructed CMB map is free from any bright point source con- 
tamination. The joint analysis also produces point source catalogues at each Planck 
frequency which are more complete and accurate than those obtained by each method 
on its own. The results are especially improved at high frequencies where infrared 
galaxies dominate the point source contribution. Although this joint technique has 
been performed on simulated Planck data, it could be easily applied to other mul- 
tifrequency CMB experiments, such as the forthcoming NASA MAP satellite or the 
recently performed Boomerang and MAXIMA experiments. 

Key words: methods: data analysis - techniques: image processing - cosmic mi- 
crowave background 



1 INTRODUCTION 

The cosmic microwave background (CMB) is one of the most 
powerful observational tools for understanding our Universe. 
Indeed, an accurate knowledge of the CMB anisotropics can 
place tight constraints on fundamental parameters such as 
the amount of matter in the Universe and its overall geome- 
try. Observations of the CMB also allow one to discriminate 
between different models of structure formation. 

Recent CMB experiments such as Boomerang (de 
Bernardis et al. 2000) and MAXIMA (Balbi et al. 2000) 
have set strong constraints on the geometry of the Universe, 
showing that it is close to spatially flat. Nevertheless, there 
remain numerous unbroken degeneracies in the full set of 
parameters that define the currently favoured inflationary 
CDM cosmological model. In order to resolve these degen- 
eracies a new generation of CMB satellite experiments is cur- 
rently in preparation, most notably the NASA MAP mission 
(Bennet et al. 1996) and the ESA Planck Surveyor (Puget et 
al. 1998, Mandolesi et al. 1998) . These experiments will pro- 
vide high-resolution all-sky maps of the CMB anisotropies 



which will allow a highly accurate estimation of a large set 
of cosmological parameters. 

An important issue for CMB satellite missions is the 
separation of foreground emission from the CMB signal. An 
accurate means of performing this separation is vital in or- 
der to make full use of the high-resolution CMB maps these 
experiments will produce. The main foregrounds to be sepa- 
rated from the CMB signal are those due to our own Galaxy 
(dust, free- free & synchrotron emission) and extragalactic 
emission due, principally, to the Sunyaev-Zel'dovich effect 
(both thermal and kinetic) and point sources. 

A preliminary application of neural-network techniques 
in this field has recently been performed with promising re- 
sults (Baccigalupi et al. 2000), but this approach is not at 
present sufficiently sophisticated to accommodate multifre- 
quency data arising from convolutions with beams of dif- 
ferent sizes and subject to different levels of instrumental 
noise. Nevertheless, more traditional techniques based on 
the maximum-entropy method (MEM) have been shown to 
provide an efficient and accurate way of performing the com- 
ponent separation (Hobson et al. 1998; H98, hereafter). 
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As one might expect, the MEM technique is partic- 
ularly successful at using multifrequency data to identify 
foreground emission from physical components whose spec- 
tral signatures are (reasonably) well-known. It is therefore 
not surprising that the most problematic foreground to re- 
move is that due to extragalactic point sources. These differ 
from the other components in the important respect that the 
spectral behaviour of point sources differs from one source 
to the next and, moreover, is notoriously difficult to predict. 
We can, however, make some headway by at least identifying 
the likely populations of point sources at different observing 
frequencies. Our knowledge of point source populations is 
increasing with new observations, but there are still great 
uncertainties. In the absence of extensive observations at 
microwave frequencies, the currently-favoured approach is 
to use the Toffolatti et al. (1998) model, which provides an 
estimate of point source populations based on existing obser- 
vations and basic physical mechanisms. This model can be 
used to generate simulated point source emission at different 
observing frequencies. From 30 to 100 GHz, the main point 
source emission is due to radio selected flat -spectrum AGNs 
(radio-loud quasars, blazars, etc.). From 300 to 900 GHz, 
infrared selected sources - starburst and late type galax- 
ies at intermediate to low redshift and high redshift ellipti- 
cals - account for most point source emission. At intermedi- 
ate frequencies, both populations contribute approximately 
equally. 

The problem of removing emission due to point sources 
from CMB observations has been addressed by several 
authors. For example, Tegmark & Oliveira-Costa (1998; 
TOC98, hereafter) suggest a straightforward harmonic fil- 
ter technique that is optimised to detect and subtract point 
sources, whereas Hobson et al. (1999; H99, hereafter) pro- 
pose treating point source emission as an additional gener- 
alised 'noise' contribution within the framework of the MEM 
approach discussed in H98. Tenorio et al. (1998) present a 
wavelet technique to subtract point sources, but the wavelet 
basis used in this work is not the optimal one for this pur- 
pose. This point is addressed in Sanz, Herranz & Martinez- 
Gonzalez (2000), where it is shown that the 'Mexican Hat' 
wavelet (MHW) is in fact optimal for detecting point sources 
under reasonable conditions, the most important assump- 
tion being that the beam is well approximated by a Gaus- 
sian. 

The application of this wavelet to realistic simulations 
has been presented in Cayon et al. (2000; COO, hereafter) and 
extended in Vielva et al. (2001; V01, hereafter) . The main 
advantage of the MHW method over the previous works is 
that the algorithm does not require any assumptions to be 
made regarding the statistical properties of the point source 
population or the underlying emission from the CMB (or 
other foreground components). 

The aim of this paper is to show how the MEM and 
MHW approaches are in fact complementary and can be 
combined to improve the accuracy of the separation of dif- 
fuse foregrounds from the CMB and increase the number 
of points sources that are identified and successfully sub- 
tracted. The technique proposed in this paper is as follows. 
First, we apply the MHW to the multifrequency simulated 
Planck maps to detect the brightest point sources at each 
observing frequency, which are subtracted from the original 
data. The MEM algorithm is then applied to these processed 



maps in order to recover the rest of the components of the 
microwave sky. These reconstructed components are then 
used as inputs to produce 'mock' data and subtracted from 
the original data maps. Since we expect our reconstructions 
to be reasonably accurate, the residuals maps obtained in 
this way would mostly contain noise and the contribution 
from point sources. Finally, we apply the MHW on these 
residuals maps in order to recover a more complete and ac- 
curate point source catalogue. 

In this paper the combined method is applied to simu- 
lated observation by the Planck satellite, but the technique 
could easily be applied to MAP data or to existing mul- 
tifrequency observations by the Boomerang or MAXIMA 
balloon experiments. The paper is organized as follows. In 
the next section we present a brief overview of the MEM 
and MHW techniques, outlining in particular the advan- 
tages and shortcomings of each approach. We then explain 
why the two approaches can be successfully combined to 
produce a more powerful joint analysis scheme. In section 3 
we summarise how the simulated Planck observations were 
performed. In section 4 we present the results of a com- 
ponent separation based on the new joint analysis method, 
and discuss the improved accuracy of the reconstructions 
of the diffuse components. In section 5 we concentrate on 
the recovery of the point sources themselves and discuss the 
construction of point source catalogues from Planck obser- 
vations. Finally, our conclusions are presented in section 6. 



2 THE MEM AND MHW TECHNIQUES 

In this section, we briefly review the MEM and MHW tech- 
niques. A complete description of the MEM component sep- 
aration algorithm can be found in H98. In addition, H99 
describes how to include point sources into the MEM for- 
malism. We therefore provide only a basic outline of the 
approach. The MHW method is introduced in COO and ex- 
tended in V01, and so again we give only a basic summary. 

2.1 The maximum-entropy method 

If we observe the microwave sky in a given direction x at 
rif different frequencies, we obtain an n/-component data 
vector that contains, for each frequency, the temperature 
fluctuations convolved with the beam in this direction plus 
instrumental noise. The ^th component of the data vector 
in the direction x may be written as 

B v (\x— x'\) F vp Sp(x') d 2 'x ' +£ v (x)+e v (x).(l) 

P =i 

In this expression we distinguish between the contributions 
from the point sources and the n c physical components for 
which it is assumed the spectral behaviour is constant and 
reasonably well-defined (over the observed patch of sky). 
The latter are collected together in a signal vector with n c 
components, such that s p (x) is the signal from the pth phys- 
ical component at some reference frequency vq. The corre- 
sponding total emission at the observing frequency v is then 
obtained by multiplying the signal vector by the rif x n c fre- 
quency response matrix F vp that includes the spectral be- 
haviour of the considered components as well as the trans- 
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mission of the uih frequency channel. This contribution is 
then convolved with the beam profile B I/ (x) of the relevant 
channel. Since the individual spectral dependencies of the 
point sources are very complicated, we cannot factorize their 
contribution in this way and so they are added into the for- 
malism as an extra 'noise' term. Thus is the contribution 
from point sources, as observed by the instrument at the 
frequency v (hence, convolved with the beam profile). Fi- 
nally, e„ is the expected level of instrumental noise in the 
z/th frequency channel and is assumed to be Gaussian and 
isotropic. 

The assumption of a spatially-invariant beam profile in 
(jjj) allows us to perform the reconstruction more effectively 
by working in Fourier space, since we may consider each k- 
mode independently (see H98). Thus, in matrix notation, at 
each mode we have 



d = Rs + | + e = Rs + C, 



(2) 



where d, £ and e are column vectors each containing Uf 
complex components and s is a column vector containing n c 
complex components. In the second equality we have com- 
bined the instrumental noise vector e and the point-source 
contribution £ into a single 'noise' vector £. The response 
matrix R has dimensions n/ x n c and its elements are given 
by R vp (k) = B v (k)F vp . The aim of any component separa- 
tion/reconstruction algorithm is to invert (^|) in some sense, 
in order to obtain an estimate s of the signal vector at each 
value of k independently. Owing to the presence of noise, and 
the fact that the response matrix R is not square and would, 
in any case, have some small eigenvalues, a direct inversion 
is not possible, and so some form of regularised inverse must 
be sought. Typical methods include singular- valued decom- 
position, Wiener filtering or the maximum-entropy method. 

The elements of the signal vector s at each Fourier mode 
may well be correlated, this correlation being described by 
the n c x n c signal covariance matrix C defined by 



C(fc) = <s(fe)s t (fc)}, 



(3) 



where the dagger denotes the Hermitian conjugate. More- 
over, if prior information is available concerning these cor- 
relations, we would wish to include it in our analysis. We 
therefore introduce the vector of 'hidden' variables h, re- 
lated to the signal vector by 



s = Lh, 



(4) 



where the n c x n c lower triangular matrix L is obtained by 
performing a Cholesky decomposition of the signal covari- 
ance matrix C = LL T . The reconstruction is then performed 
entirely in terms of h and the corresponding reconstructed 
signal vector is subsequently found using (Q). 

Using Bayes' theorem, we choose our estimator h of 
the hidden vector to be that which maximises the posterior 
probability given by 



Pr(h|d) tx Pr(d|h)Pr(h) 



(•») 



where Pr(d|h) is the likelihood of obtaining the data given 
a particular hidden vector and Pr(h) is the prior probability 
that codifies our expectations about the hidden vector before 
acquiring any data. 



As explained in H99, the form of the likelihood function 
in (|E|) is given by 

Pr(d|h) tx exp(-C t N _1 C) 

oc exp[-(d-RLh) t N" 1 (d-RLh)] (6) 

where in the last line we have used (Q) . The noise covariance 
matrix N has dimensions nf x nf and at any given fc-mode 
is given by 



N(fc) = (C(*)C t (*)>- 



(7) 



Therefore, at a given Fourier mode, the ^th diagonal element 
of N contains the ensembled- averaged power spectrum at 
that mode of the instrumental noise plus the point source 
contribution to the vth frequency channel. The off-diagonal 
terms give the cross-correlations between different channels; 
if the noise is uncorrelated between channels, only the point 
sources contribute to the off-diagonal elements. 

For the prior Pr(h) in (ph, we assume the entropic form 



Pr(h) oc exp[aS'(h,m) 



(8) 



where S(h, m) is the cross entropy of the complex vectors 
h and m, where m is a model vector to which h defaults in 
absence of data. The form of the cross entropy for complex 
images and the Bayesian method for fixing the regularising 
parameter a are discussed in H98. We note that, in the ab- 
sence of non-Gaussian signals, the entropic prior (Q) tends 
to the Gaussian prior implicitly assumed by Wiener filter 
separation algorithms, and so in this case the two methods 
coincide. 

The argument of the exponential in the likelihood func- 
tion (^) may be identified as (minus) the standard \ 2 misfit 
statistic, so we may write Pr(d|h) oc exp[— x 2 (h)]. Substitut- 
ing this expression, together with that for the prior proba- 
bility given in (^|), into Bayes' theorem, we find that max- 
imising the posterior probability Pr(h|d) with respect to h 
is equivalent to minimising the function 

*(h) = X 2 (h)-«S(h,m). 

This minimisation can be performed using a variable metric 
minimiser (Press et al. 1994) and requires only a few minutes 
of CPU time on a Sparc Ultra workstation. 



2.2 The mexican hat wavelet method 

The MHW technique presented by COO and V01 for identi- 
fying and subtracting point sources operates on individual 
data maps. Let us consider the two-dimensional data map 
d v (x) at the frequency v. If the map contains point sources 
at positions bi with fluxes or amplitudes Ai, together with 
contributions from other physical components and instru- 
mental noise, then the data map is given by 



d v {x) = £ v (x) + n v {x) = AiB v (x — bi) + n v (x 



(9) 



where B v {x) is the beam at the observing frequency v, and 
in this case the generalised 'noise' n v (x) is defined as all 
contributions to the data map aside from the point sources. 

For the ith point source, we may define a 'detection 
level' 



D = 



At/Q 



(10) 
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where Q. is the area under the beam and a n is the dispersion 
of the generalised noise field n v (x). In general, the detec- 
tion level D will be much less than unity for all but the 
few brightest sources. This is the usual problem one faces 
when attempting to identify point sources directly in the 
data map. 

As explained in COO, instead of attempting the iden- 
tification in real space, one can achieve better results by 
first transforming to wavelet space. For a two-dimensional 
data map d v (x), we define the continuous isotropic wavelet 
transform by 



Wd{R, b) 



d v (x), 



(11) 



where w(R, b) is the wavelet coefficient associated with the 
scale R at the point b (where the wavelet is centred). The 
function VKM) is t ne mother wavelet, which is assumed to 
be isotropic and satisfies the conditions 



fd 2 xip(x) = 

Jd 2 x\i}(x)\ 2 = 1 

CV = (2tt) 2 J™ dkk- 1 ]^)] 2 < 



(compensation) , 
(normalisation) , 
(admissibility), 



where tp denotes the Fourier transform of x = \x\ and 
k — \k\. The wavelet coefficients given by (hit) characterise 
the contribution from structure on a scale R to the value of 
the map at the position b. 

By analogy with (|Io|), in wavelet space we define the 
detection level for the ith point source (as a function of 
scale R) by 



D W (R) = 



w l (R, bj) 
a Wn (R) 



(12) 



where w^(R, bi) is the wavelet coefficient of the field £„(x) 
at the location of the ith point source, and a 2 Un (R) is the 
dispersion of the wavelet coefficients w„(R, b) of the gener- 
alised noise field n v (x). It is straightforward to show that 



al n {R) = 2'K 



dk kP{k)\il){Rk)Y 



(13) 



where P(k) is the power spectrum of n v (x). The integral 
limits, k m i n and k max , correspond to the maximum and 
minimum scales of the sky patch analysed, i.e. the patch 
and pixel scales respectively. 

The detection level D w (R) in wavelet space will have a 
maximum value at some scale R — Rq. This scale is prac- 
tically the same for all point sources, and may be found by 
solving dD w (Ro)/dR — 0. In general, the optimal scale is of 
the order of the beam dispersion a a (see V01, section 3, for 
a discussion about how the noise and the coherence scale of 
the background determine the optimal scale). In order that 
the wavelet coefficients are optimally sensitive to the pres- 
ence of the point source, we must make the value of D w (Ro) 
as large as possible. This is achieved through an appropri- 
ate choice both of the mother wavelet ip(x) and the optimal 
scale. If the beam profile is Gaussian and the power spec- 
trum of the generalised noise field is scale-free, Sanz et al. 
(2000) show that for a wide range of spectral indices of the 
power spectrum the Mexican Hat wavelet is optimal. The 
two-dimensional MHW is given by 



from which we find 



W((R, bi) = 2\/2^ 



A (l+(i?/<7 a )2)2 : 



(14) 



(15) 



where At is the amplitude of the ith point source, Q is the 
area under the Gaussian beam and a a is the beam disper- 
sion. In (|l5|) it is assumed that any overlap of the (convolved) 
point sources is negligible. That is a good approximation for 
the brightest point sources, the ones that the MHW is able 
to detect. In fact, the number of point sources detected at 
each frequency (see Table |^) corresponds only to a small 
percentage of the number of resolution elements (see Ta- 
ble contained at each Planck frequency channel (~ 6% at 
30GHz, the most unfavourable case). Therefore, the proba- 
bility of finding two or more bright point sources inside the 
same resolution element is very low. The advantage of iden- 
tifying point sources in wavelet space rather than real space 
may then be characterised by the amplification factor 



A = 



D w (Ro) 
D 



2V2-K 



(Ro/c 



(1 + {R /(J a YY a Wn (Ro)' 



(16) 



In practice, it is clear that we do not have access to 
the wavelet coefficients of the fields £v(x) and n v (x) sepa- 
rately, but only to the wavelet coefficients of the total data 
map d v (x) = £v(x) + n v (x). Nevertheless, if the detection 
level D W (R) for the ith point source is reasonably large, we 
would expect w^(R, bi) « Wd{R, bi). Also, if we assume that 
the power spectrum of the point source emission is negligi- 
ble as compared to that of the generalised noise field, then 
o~w n (R) ~ cr Wd (R). Thus our algorithm for detecting point 
sources is as follows. Using the above approximations, we 
first calculate the optimal scale Ro- We then calculate the 
wavelet transform Wd(Ro, b) of the data map at the opti- 
mal scale. The wavelet coefficients are then analysed to find 
sets of connected pixels above a certain threshold a Wd (Ro). 
The maxima of these spots are taken to correspond to the 
locations of the point sources. 

For every point source detected in the above way, we 
then go on to estimate its flux. This is achieved by perform- 
ing a multiscale fit as follows. For each point source location 
bi, the wavelet transform Wd{R, bi) is calculated at a number 
of scales R and compared with the theoretical curve (|l5[). 
This comparison is performed by calculating the standard 
misfit statistic 



X 



„( cx p) 



(theo)] 



,,( ex p) 



, (theo)"] 



„( gx p) 



(17) 
(Rk,b t ) 



where the fcth element of the vector ™ (oxp) is 
(and similarly for the vector of theoretical wavelet coeffi- 
cients). The matrix V is the empirical covariance matrix of 
the wavelet coefficients on different scales, which is given by 



where the average is over position b. 



(18) 



2.3 MEM and MHW joint analysis 

In H99 the MEM technique is shown to be effective at per- 
forming a full component separation in the presence of point 
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sources. In particular, the reconstructed maps of the sep- 
arate diffuse components contain far less point-source con- 
tamination than the input data maps. Moreover, by compar- 
ing the true data maps with simulated data maps produced 
from the separated components, it is possible to obtain point 
source catalogues at each observing frequency. Nevertheless, 
the MEM approach does have its limitations. Since the point 
sources are modelled as an additional generalised noise com- 
ponent, it is not surprising that MEM performs well in iden- 
tifying and removing the large number of point sources with 
low to intermediate fluxes. However, it is rather poorer at 
removing the contributions from the brightest point sources. 
These tend to remain in the reconstructed maps of the sep- 
arate diffuse components, although with much reduced am- 
plitudes. 

The MHW technique, on the other hand, performs best 
when identifying and removing the brightest point sources. 
Indeed, in detecting bright sources the MHW technique 
generally out-performs other techniques such as SExtractor 
(Bertin & Arnouts, 1996) and standard harmonic filtering 
(TOC98). Moreover, the amplitudes of the bright sources 
are also accurately estimated. For weaker sources, however, 
the MHW technique performs more poorly by either inac- 
curately estimating the flux or failing to detect the source 
altogether. 

The strengths and weaknesses of the MEM and MHW 
approaches clearly indicate that they are complementary 
techniques, and that a combined approach might lead to im- 
proved results as compared to using each method indepen- 
dently. We thus propose the following method for analysing 
multifrequency observations of the CMB that contain point 
source contamination. First, the data map at each observing 
frequency is analysed separately using the MHW in order to 
detect and remove as many bright point sources as possible 
and obtain accurate estimates of their fluxes. The processed 
data maps are then taken as the inputs to the generalised 
MEM approach discussed in H99. As we will demonstrate in 
section Q the MEM analysis of these processed maps leads 
to more accurate reconstructed maps of the separate diffuse 
components. This leads in turn to more accurate residual 
maps between the true input data with the data simulated 
from the reconstructions. These residual maps are then anal- 
ysed with the MHW in order both to refine the original 
estimates of the fluxes of the bright point sources and to de- 
tect fainter sources. Thus, the joint analysis not only gives 
a more complete point source catalogue, but also improves 
the quality of the reconstructed maps of the CMB and other 
foreground components. 



3 SIMULATED OBSERVATIONS 

As mentioned in the introduction, the joint analysis tech- 
nique can be applied to any multifrequency observations of 
the CMB that may contain point source emission. Thus, for 
example, the method could straightforwardly be applied to 
the existing Boomerang or MAXIMA data-sets, or to ob- 
servations by the forthcoming NASA MAP satellite. Never- 
theless, in order to test the capabilities of the joint analysis 
method to the fullest extent, in this paper we apply it to sim- 
ulated observations by the proposed ESA Planck Surveyor 



Table 1. Basic observational parameters of the 10 frequency 
channels of the Planck Surveyor satellite. Column two lists the 
fractional banwidths. The FWHM in column three assumes a 
Gaussian beam. In column four the instrumental noise level is 
AT (/iK) per resolution element for 12 months of observation. 



Frequency 
(GHz) 


Fractional bandwidth 
(Av/u) 


FWHM 
(arcmin) 


G noise 

GuK) 


30 


0.20 


33.0 


4.4 


44 


0.20 


23.0 


6.5 


70 


0.20 


14.0 


9.8 


100 (LFI) 


0.20 


10.0 


11.7 


100 (HFI) 


0.25 


10.7 


4.6 


143 


0.25 


8.0 


5.5 


217 


0.25 


5.5 


11.7 


353 


0.25 


5.0 


39.3 


5 15 


0.25 


5.0 


400.7 


857 


0.25 


5.0 


18182 



satellite. The basic observational parameters of the Planck 
mission instruments (HFI and LFI) are listed in Table |l| 

The simulated data are similar to those used in H99. 
They include contributions from the primordial CMB, the 
thermal and kinetic Sunyaev-Zel'dovich effects, extragalac- 
tic point sources and Galactic thermal dust, free-free and 
synchrotron emission. Aside from the point sources, we as- 
sume that the emission of each physical component can 
be factorised into a spatial template at 300 GHz with a 
known frequency dependence. In Figure |^ we have plotted 
the six component templates at 300 GHz. Each map covers 
a 12.8 x 12.8 deg 2 patch of sky and has been convolved with 
a Gaussian beam with FWHM 5 arcmin (i.e. the highest 
Planck resolution). 

The CMB map is a Gaussian realisation of a spatially- 
flat inflationary/CDM model with Q. m = 0.3 and Q.\ = 0.7, 
for which the Ci coefficients were generated using CMB- 
FAST (Seljak & Zaldarriaga 1996). The thermal Sunyaev- 
Zel'dovich map is taken from the simulations of Diego et 
al. (2000) and assume the same cosmological model as that 
used for the CMB. The kinetic SZ field is produced by as- 
suming the line-of-sight cluster velocities are drawn from a 
Gaussian distribution with zero mean and rms 400 km/s. 
The extragalactic point source simulations adopt the model 
of Toffolatti et al. (1998) and also assume the same cosmo- 
logical model. 

The Galactic thermal dust emission is created using 
the template of Finkbeiner, Davis & Schlegel (1999). The 
frequency dependence of the dust emission is assumed to 
follow a grey-body function characterised by a dust temper- 
ature of 18 K and an emissivity /3 = 2. The distribution of 
Galactic free-free emission is poorly known. Current experi- 
ments such as the H-a Sky Survey^] and the WHAM project^] 
should soon provide maps of H a emission that could be used 
as templates. For the time being, however, we create a free- 
free template that is correlated with the dust emission in the 
manner proposed by Bouchet, Gispert & Puget (1996). The 
frequency dependence of the free-free emission is assumed 

* http:/ /www. swarthmore.edu/Home/News/Astronomy/ 
t http:/ /www. astro. wisc.edu/wham/ 
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Figure 1. The rms thermodynamic temperature fluctuations at 
the Planck observing frequencies due to each physical component, 
after convolution with the appropriate beam and asumming a 
sampling rate of FWHM/2.4. The rms instrumental noise per 
pixel at each frequency is also plotted. 



to vary as /„ oc u~ 016 , and is normalised to give an rms 
temperature fluctuation of 6.2fiK at 53 GHz. Finally, the 
synchrotron spatial template has been produced using the 
all sky map of Fosalba & Giardino[|]. This map is an extrap- 
olation of the 408 MHz radio map of Haslam et al. (1982), 
from the original 1 deg resolution to a resolution of about 
5 arcmin. The additional small-scale structure is assumed 
to have a power-law power spectrum with an exponent of 
—3. We have continued this extrapolation to 1.5 arcmin fol- 
lowing the same power law. The frequency dependence is 
assumed to be I v oc is~ 9 and is normalized to the Haslam 
408 MHz map. 

To simulate the observed data in a given Planck fre- 
quency channel, each of the physical components discussed 
above is first projected to the relevant frequency and the 
contributions are summed. The predicted point source emis- 
sion for the frequency is then added, and the resulting total 
sky emission is convolved with a Gaussian beam of the ap- 
propriate FWHM. Finally, independent pixel noise is added, 
with corresponding rms from Table In Figure [l] we give the 
rms thermodynamic temperature fluctuations in the data at 
each Planck observing frequency due to each physical com- 
ponent and the instrumental noise. In Figure H we plot the 
simulated Planck observations in each frequency channel, 
all the components are included: CMB, dust, free-free, syn- 
chrotron, kinetic and thermal SZ effects and point source 
emission as well as instrumental noise. 



4 FOREGROUND SEPARATION 



We have applied the method outlined in section 2.3 to the 
simulated Planck data described above. We have assumed 
knowledge of the azimuthally averaged power spectra of 
the six input components in Fig. |^, together with the az- 
imuthally averaged cross power spectra between them (see 
H98 for more details). Using the model of Toffolatti et al. 

"t ftp:/ /astro. estec.esa.nl/pub/synchrotron 



Table 2. The rms in /iK of the reconstruction residuals smoothed 
with a 5 arcmin FWHM Gaussian beam with and without the 
initial subtraction of bright point sources using the MHW. Full 
power spectrum information has been assumed. For comparison 
the rms of the input maps shown in Fig. [| are also given. Results 
are given for the reference frequency of 300 GHz. 



Component 


input 


error 


error 




rms 


(with MHW) 


(without MHW) 


CMB 


112.3 


7.68 


8.62 


Kinetic SZ 


0.69 


0.70 


0.70 


Thermal SZ 


5.37 


4.64 


4.66 


Dust 


55.8 


2.68 


3.39 


Free- Free 


0.66 


0.22 


0.24 


Synchrotron 


0.32 


0.11 


0.12 



(1998), we have also introduced the power spectrum of the 
point sources at each frequency channel, including cross 
power spectra between channels, and account for this con- 
taminant as an extra noise term (see H99 for more details). 
However, the recovery of the main components and point 
sources does not depend critically on this assumption, as 
will be discussed later. 

The resulting reconstructions of the physical compo- 
nents at a reference frequency of 300 GHz are shown in 
Fig. ^. The maps have been plotted using the same grey- 
scale as in Fig. ^ to allow a straightforward comparison. In 
Fig. [5], we plot the residuals for each component, obtained 
by subtracting the input maps from the reconstructions. 

We can see that the main input components have been 
faithfully recovered with no obvious visible contamination 
from point sources. This is because the MHW subtraction 
algorithm is efficient at removing the brightest point sources, 
whereas MEM has greatly reduced the contamination due 
to fainter sources. We give the rms reconstruction errors for 
each component in Table ^| 

In particular, we note that the CMB has been recovered 
very accurately, although the residuals map does show some 
weak contamination due to low-amplitude point sources. In- 
deed, the rms reconstruction error for this component is ~ 
7.7 pK, which corresponds to an accuracy of ~ 6.8 per cent 
as compared to the rms of the input CMB map (see Table |J). 
Even more impressive is the reconstruction of the dust map. 
None of the numerous point sources present in the highest 
frequency channel maps are visible in the reconstruction. 
This is also confirmed by inspecting the residual map. The 
main features of the free-free emission are also recovered, 
mostly due to its high correlation with the dust. Again, the 
reconstruction shows no evidence of point source contamina- 
tion. For the synchrotron component, the recovered emission 
is basically a lower resolution image of the input map. This 
is expected since only the lowest frequency channels provide 
useful information about this component, and these channels 
also have the lowest angular resolutions. Although the recon- 
structed synchrotron map is mostly free of point sources, 
some residual contamination remains. This contamination 
corresponds to a few medium amplitude point sources that 
are present in the lowest frequency channels, although they 
are not clearly visible in the data. These sources are too 
weak to be detected using the MHW algorithm but at the 
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Figure 2. The 12.8 X 12.8 deg 2 realisations of the six input components used to produce the simulated Planck data. The different 
panels correspond to (a) CMB, (b) kinetic SZ effect, (c) thermal SZ effect, (d) Galactic dust, (e) Galactic free-free and (f) Galactic 
synchrotron emission. Each component is plotted at 300 GHz and has been convolved with a Gaussian beam of FWHM 5 arcmin (the 
highest resolution expected for the Planck satellite). The map units are equivalent thermodynamic temperature in fiK. 



same time they are not well characterised by the generalised 
noise approach assumed in MEM. 

As pointed out in H99, one must be careful when com- 
paring the amplitude of the residual point sources still con- 
taminating the reconstructions with the corresponding am- 
plitudes of the point sources in the data maps. The recon- 
structions are calculated at a reference frequency of 300 GHz 
and those sources remaining in the residuals maps are pro- 
jected in frequency according to the spectral dependence 
of the component they contaminate. In addition, we have 



to take into account the different resolution of the Planck 
frequency channels. For example, the contaminating point 
source in the middle right-hand side of the synchrotron resid- 
uals map has an amplitude of ~ 0.15/iK after convolution 
with a Gaussian beam of FWHM 33 arcmin (the resolu- 
tion of the 30 GHz Planck frequency channel). Following 
the spectral dependence of the synchrotron component, this 
projects to 17.5^iK at 30GHz. This value should be com- 
pared with the amplitude of the point source at the same 
frequency, which is around 150piK. Therefore, MEM has 
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Figure 3. The 12. 8x 12.8 deg 2 data maps observed at each of the ten Planck frequency channels listed in Table|lJ The panels correspond to 
the frequencies (a) 30 GHz, (b) 44 GHz, (c) 70 GHz, (d) 100 GHz-lfi, (e) 100 GHz-hfi, (f) 143 GHz, (g) 217 GHz, (h) 353 GHz, (i) 545 GHz 
and (j) 857 GHz. The units are equivalent thermodynamic temperature in /jK. 
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Figure 4. The reconstructed maps for each of the physical components: (a) CMB, (b) kinetic SZ effect, (c) thermal SZ effect, (d) 
Galactic dust, (e) Galactic free-free and (f) Galactic synchrotron emission. Point sources have been subtracted from the data maps using 
the mexican hat algorithm before applying MEM. Each component is plotted at 300 GHz and has been convolved with a Gaussian beam 
of FWHM 5 arcmin. The map units are equivalent thermodynamic temperature in p,K. 



succeeded in reducing the contamination due to this point 
source by almost a factor of ten. 

The recovery of the thermal SZ effect is quite good. 
Most of the bright clusters have been reproduced whereas 
only a few point sources has been misidentified as clusters. 
At the reference frequency of 300GHz, these misidentified 
point sources appear mostly as negative features. Finally, as 
expected, the reconstruction of the kinetic SZ is quite poor 
and one detects only a few clusters whose corresponding 
thermal SZ effect is large. 



We have also calculated the power spectrum of the re- 
constructed component maps and found that the accuracy 
is very similar to that found in H99, so we do not plot them 
again here. The effect of first applying the MHW to the data 
maps before the MEM analysis is not so obvious when con- 
sidering the power spectra of the reconstructions, since only 
a small percentage of pixels are affected by residual point 
sources and this has little effect in the recovered spectrum. 
Nevertheless, the removal of the point source contamination 
is vital if one wishes to probe the Gaussian character of the 
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Figure 5. The reconstruction residuals obtained from subtracting the input maps of Fig. g from the reconstructed maps of Fig. ^. 
The panels correspond to (a) CMB; (b) kinetic SZ effect; (c) thermal SZ effect; (d) Galactic dust; (e) Galactic free-free; (f) Galactic 
syncrotron emission. 



CMB, as well as to study properties of the other foreground 
components. 

To understand better the effect of performing the MHW 
analysis on the data maps prior to the MEM algorithm, we 
have also carried out a component separation for the case 
where the MHW step is not performed; this corresponds 
to the method in H99. In Fig. ^, we show the difference 
between the reconstructions obtained using the combined 
MHW and MEM technique and those obtained using MEM 
alone. Thus these maps display the point sources that have 
been successfully removed by the MHW, which would be 



otherwise present in the reconstructions. By comparing with 
the data, we can also see how the point sources present in 
the different frequency channels would affect a given com- 
ponent if not carefully subtracted. Particularly impressive 
is the removal from the dust and free-free reconstructions 
of the large number of point sources that were present in 
high-frequency data channels. For the CMB, the MHW has 
subtracted a few very bright point sources, which dominated 
the contribution of this contaminant in the intermediate fre- 
quency channels. We also note that the MHW has removed 
a few from the synchrotron reconstruction that were present 
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Figure 6. The residuals obtained by subtracting the reconstructed maps shown in Fig. [jj from the reconstructions achieved when the 
MHW is not used to perform an initial point source removal. Both reconstructions were smoothed with a Gaussian beam of FWHM 5 
arcmin and correspond to a frequency of 300 GHz. The different panels are: (a) CMB, (b) kinetic SZ effect, (c) thermal SZ effect, (d) 
Galactic dust, (e) Galactic free-free and (f) Galactic synchrotron emission. The map units are equivalent thermodynamic temperature 
in /iK. 



in the lowest channels. The reconstructions of the thermal 
and kinetic SZ effects have also been improved since a lower 
number of point sources have been misidentified as clusters; 
these sources were detected by the MHW mainly in the high- 
est channels of the LFI. The rms reconstructions errors when 
MEM is used without a previous subtraction of point sources 
by the MHW are also given in Table ^. 

Finally, we can also study how our reconstructions are 
affected if we do not assume full power spectrum informa- 



tion. Thus, we have repeated our joint analysis of the simu- 
lated Planck observations for the case where we assume that 
£ 2 Ce is constant for each component out to the highest mea- 
sured Fourier mode. The level of the fiat power spectrum for 
each component is, however, chosen so that the total power 
in each component is approximately that observed in the in- 
put maps in Fig. ^. Furthermore, no information about the 
cross power spectra between different components is given. 
Regarding the point sources, the true azimuthally averaged 
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power spectrum is again assumed to account for their con- 
tribution as an extra noise term but cross-correlations be- 
tween different frequency channels are ignored. The quality 
of the reconstructions of the main components is actually 
very similar to the case when full power spectrum informa- 
tion is given. In particular, the accuracy of the CMB re- 
constructed map is only slightly worse with a rms error of 
8.2/iK as compared to 7.7/iK in the former case. Moreover, 
the reconstruction is again free from obvious contamination 
due to point sources. Similarly, the dust component has been 
faithfully recovered with a rms error of 3.0pK versus 2.7/iK 
when full power spectrum is assumed. The main features 
of the synchrotron and thermal SZ effect are also recovered 
although the reconstructions are poorer. In particular, the 
contamination due to point sources is slightly increased in 
the thermal SZ reconstructed map. Finally, the weakest com- 
ponents, free- free emission and kinetic SZ effect, are lost in 
this case and the reconstructions have simply defaulted to 
zero in the absence of any useful information. 



5 RECOVERY OF POINT SOURCES 

The previous section focused on the recovery of the six phys- 
ical component maps shown in Fig. |j| However, a major aim 
of the Planck mission is also to compile point source cata- 
logues at each of its observing frequencies. In this section, 
we compare the catalogues obtained using MEM alone (i.e., 
without a previous subtraction of bright sources detected 
by the MHW), MHW alone and the joint analysis method 
M&M. 

The MHW catalogue (MHWc) is produced in the man- 
ner explained in V01, through the so called 50% error cri- 
terion (see the previous work for more details in detection 
criteria). In short, each of the data maps of Fig. |^ is inde- 
pendently analysed with the MHW. Those coefficients above 
2a Wd at the optimal wavelet scale are identified as point 
source candidates. A multiscale fit is then performed in or- 
der to estimate their amplitude and those wavelet coeffi- 
cients with a non-acceptable fit are discarded. We then look 
for the flux above which at least 95 per cent of the the de- 
tected point sources have a relative error ^ 50%. This gives 
our estimation of the flux limit. Thus, the number of de- 
tections at a given channel is given by those point sources 
with an estimated amplitude above the flux limit. In prac- 
tice, we also use multifrequency information to include those 
point sources that, having an error larger than 50% or an 
insufficiently good fit, have been detected in an adjacent 
channel (although in most channels this only accounts for a 
very small fraction of the detected point sources). The error 
is defined as: E — \A — A e \/A, where A is the flux of the 
simulated source and A e that of the estimated one. 

Although the MHW catalogue (hereafter MHWc) is ob- 
tained in the same way as that of V01, the results here differ 
slightly. This apparent discrepancy is due to the different 
sampling rates that have been considered in each case. V01 
assumes pixel sizes of 1.5, 3 and 6 arcminutes for the differ- 
ent Planck channels, whereas in this work the pixel size is 
given by a fixed sampling of 2.4, following by a regridding 
to pixels of size 1.5 arcminutes. Therefore, the number of 
detected point sources and fluxes of V01 and this paper are 
not directly comparable. 



Nevertheless, not all the point sources in the MHWc 
are subtracted from the original maps prior to performing 
the MEM analysis. In theory, giving as much information as 
available should improve the MEM results. However, when 
using the 50% error criterion a significant number of point 
sources, especially at the 545 and 857 GHz channels, are es- 
timated with a large error. Therefore, if this information is 
given to MEM (by subtracting these sources from the origi- 
nal maps) , we are misleading the MEM algorithm. To avoid 
this unwanted effect, the point sources subtracted from the 
maps should be those with the lowest errors in the ampli- 
tude estimation. We need a more robust criterion than that 
of the 50% error. Instead we adopt the so called 5a Wd cri- 
terion which is also explained in VOL Briefly, we consider 
that a point source has been detected if the position of its 
maximum is above 5a Wd and its multiscale fit is acceptable. 

The MEM catalogue (MEMc) and the M&M one 
(M&Mc) are obtained using the method outlined in H99. 
First, the reconstructed maps of Fig. ^ are used as inputs 
to produce 'mock' data. We follow the same procedure as 
that used to obtain the data of Fig. but, of course, we do 
not add instrumental noise or the point sources. These mock 
data are then subtracted from the true data (which contain 
the full point source contribution). Since the reconstructions 
of the six main components are reasonably accurate and also 
contain very few point sources, we obtain a data residuals 
map at each Planck frequency that consists mainly of the 
point sources and instrumental noise. Each of these residuals 
maps are then independently analysed in order to produce a 
point souce catalogue at each observing frequency. We point 
out, however, that the residuals maps produced here differ 
from those in H99. In order to concentrate on the effect of 
emission from other physical components on the point source 
recovery, in H99 the instrumental noise was neglected when 
making the residuals maps. Here the instrumental noise is 
included to obtain a more realistic estimate of the number 
of points sources recoverable from real Planck data. 

Another difference with H99 is the process by which 
the point source catalogue is produced from the residuals 
maps. In H99, the SExtractor package (Bertin & Arnouts, 
1996) is used to detect and estimate the amplitude of point 
sources. The SExtractor package begins by fitting and sub- 
tracting an unresolved background component, before iden- 
tifying any point sources, and can lead to ambiguities in 
assigning a flux detection limit. Therefore, in the present 
paper, the residuals maps are instead analysed using the 
MHW, since this wavelet filter is optimal for this purpose 
(Sanz et al. 2000). At this point, however, it is worth point- 
ing out some subtleties associated with applying the MHW 
in these circumstances. In particular, for the MEMc and 
M&Mc, we apply the 50% error criterion into the residual 
maps, in order to compare with the MHWc. Clearly, this 
choice determines empirically the flux limit of the catalogue 
(achieved by the 50% error criterion) and will depend on 
the assumed point source population model, but we expect 
that most models lead to similar results. In fact, in V01 it 
is shown that for the Guiderdoni et al. (1998) E model, the 
flux limits achieved are very similar. 

In Table ^| we give, for each catalogue, the number of 
point sources detected, the minimum flux achieved and the 
average amplitude error in each Planck frequency channel. 

In the two highest frequency channels, the joint anal- 
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Table 3. The point source catalogues obtained using the MHW alone (MHWc), MEM alone (MEMc) and the joint analysis method 
(M&Mc). For each Planck observing frequency, we list the number of sources detected, the flux limit of the catalogue and the average 
of the absolute value of the percentage error for the amplitude estimation. The numbers in brackets in the second column indicate 
the number of point sources that are detected by the 5(T Wd criterion (see text for details). 







MHWc 






TV yTTTlA It 

MEMc 






A T AT 

M&Mc 




Frequency 


Number of 


Min Flux 


E 


Number of 


Min Flux 


E 


Number of 


Min Flux 


E 


(GHz) 


detections 


(Jy) 


(96) 


detections 


(Jy) 


(%) 


detections 


(Jy) 


(%) 


30 


15 (4) 


0.13 


14.0 


32 


0.09 


14.9 


30 


0.09 


14.0 


44 


11 (3) 


0.27 


11.2 


28 


0.11 


14.2 


28 


0.11 


13.2 


70 


7(5) 


0.30 


7.9 


38 


0.08 


11.8 


35 


0.09 


11.4 


100 (LFI) 


12 (3) 


0.23 


11.8 


37 


0.08 


16.4 


44 


0.06 


18.4 


100 (HFI) 


17 (7) 


0.14 


14.1 


72 


0.04 


17.3 


74 


0.03 


17.4 


143 


11 (4) 


0.16 


18.0 









5 


0.32 


12.4 


217 


15 (5) 


0.08 


14.4 









5 


0.25 


15.4 


353 


16 (10) 


0.15 


14.3 


38 


0.08 


28.7 


37 


0.08 


34.9 


545 


37 (29) 


0.29 


14.1 


89 


0.18 


26.4 


121 


0.15 


17.2 


857 


306 (86) 


0.30 


16.5 


458 


0.23 


20.6 


492 


0.19 


19.5 



ysis clearly outperforms the results obtained with each of 
the methods separately. This is due to the complementary 
nature of the two approaches, so that bright sources are de- 
tected by MHW and intermediate flux sources are identified 
by MEM. If MEM alone is used, many of the brightest point 
sources remain in the MEM reconstructions since they are 
not well characterised by a generalised noise. Therefore they 
are either not detected in the data residuals maps or the 
error in the estimated amplitude is significantly large. How- 
ever, in the joint analysis these sources are easily detected 
and their fluxes accurately estimated. 

Regarding the 353 GHz channel, the number of point 
sources detected with M&M is comparable to the best of the 
individual methods, i.e., when using MEM on its own. This 
is due to the fact that the main contaminant of the residuals 
maps is the high level of noise of the 353 GHz channel, which 
is equally present in both the residuals obtained with MEM 
and with M&M. In fact, many of the point sources have been 
detected (with a large error) thanks to the multifrequency 
information. 

The low number of point sources detected at interme- 
diate frequencies (217 and 143 GHz) in the M&Mc in com- 
paration with the MHWc, can be explained as a combina- 
tion of factors. First of all, due to the high level of noise in 
these channels relative to the point source emission, only a 
few point sources are detected with the MHW in the first 
step of our analysis, through the 5a Wd criterion. Therefore, 
when applying MEM, part of the emission of the undetected 
point sources is left in the reconstructed components, being 
mainly misidentified with CMB, which is the dominant com- 
ponent at those frequencies. This has the effect of lowering 
the level of the point sources in the data residuals, which 
together with the high level of noise, leads to a low number 
of recovered point sources. 

In the low-frequency channels, the number of point 
sources detected by the joint analysis is similar to that by 
MEM alone (the fact that M&M does not work better than 
MEM alone for all the channels may be understood as a 
statistical fluctuation; we expect that with an all sky point 
source simulation, the results of the combined method will 
be better than those of MEM in every case). In this case, the 



MHW contribution is to improve the amplitude estimation 
of a few point sources, which leads to a lower average ampli- 
tude error. In these frequency channels MHW alone detects 
only the few brightest point sources. This is because these 
channels are dominated by the CMB and the beam FWHM 
is so large that the CMB and point sources have a similar 
characteristic scale. 

Regarding the estimation of point source amplitudes, 
the joint analysis also performs better than each method 
independently. Although the average error in the amplitude 
estimation can be larger in the M&Mc than in the MHWc 
due to the detection of a larger number of faint point sources, 
those point sources present in all three catalogues are, on 
average, better estimated with the combined analysis. 

In Figure |^ we plot the amplitude estimation errors for 
the MHWc, MEMc and M&Mc versus the true flux (in Jy) 
for two representative channels: 44 and 545 GHz. We can see 
that there is a clear bias in the estimation of the amplitude 
of the brightest point sources in MEMc since they remain in 
the reconstructions and are therefore underestimated (corre- 
sponding to positive errors in Fig. . This problem is solved 
when combining MEM and the MHW. It is also obvious 
that a larger number of point sources and fainter fluxes are 
achieved in the combined analysis with respect to the MHW 
on its own. 

In Figure ^ we have plotted the sources in M&Mc, for 
the same two channels, together with the input point sources 
maps. We see that the main features are very well recovered. 
At 44 and 545 GHz, the flux limit is comparable to the level 
of instrumental noise (see next Section). Thus, to increase 
still further the number of point sources detected and reach 
fainter fluxes, one would need to denoise the residuals maps; 
this is discussed in the next section. 

We have also investigated the effect of reducing the 
power spectrum information given to MEM. We find that 
even in the extreme case when a flat power spectrum is 
assumed for the different components, the results are not 
significantly different for the M&Mc, but the quality of the 
MEMc is somewhat reduced, especially in the high frequency 
channels. In particular, the amplitude estimation errors are 
higher and the catalogue flux limit increases slightly. 
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Error in the amplitude estimation at 44 GHz 



100 
80 h 
60 
40 
20 


-20 
-40 
-60 
-80 
-100 



v v 

□ft 



V 



^ 5? 

* 



* 

* 

□ 
V 



□ 
V 

* 



* 



* MHWc 
□ MEMc 
V M&Mc 



0.1 



0.25 



0.5 



0.75 



Flux (Jy) 



100 
80 
60 
40 
20 


-20 
-40 
-60 
-80 
-100 



Error in the amplitude estimation at 545 GHz 
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Figure 7. Errors in the amplitude estimation for the MHW, MEM and joint analysis catalogues (MHWc, MEMc and M&Mc, respectively) 
as a function of the flux. We plot two Planck frequencies: 44 GHz (top) and 545 GHz (bottom). 



6 DISCUSSION AND CONCLUSIONS 

The MEM (H98, H99) and the MHW (COO, V01) techniques 
have complementary characteristics when recovering the mi- 
crowave sky. On the one hand, the MEM technique is a 
powerful tool for using multifrequency data to separate the 
cosmological signal from foreground emission whose spectral 
behaviours are (reasonably) well-known. The most problem- 
atic foreground to remove is that due to point sources. On 
the other hand, the MHW has shown to be a robust and self- 
consistent method to detect and subtract this point source 
emission from microwave maps. The aim of this paper has 
been to show how the performance of a combined (MEM and 



MHW) analysis can improve the recovery of the components 
(CMB, Sunyaev-Zel'dovich, extragalactic point sources and 
Galactic emission) of simulated microwave maps. In order 
to test this analysis, we have applied it to simulated ESA 
Planck satellite observations. However, the technique could 
straightforwardly be applied to other CMB experiments (e.g. 
NASA MAP satellite, Boomerang and MAXIMA). 

The proposed method to analyse these data is as fol- 
lows. First, we apply the MHW at each observing frequency 
in order to remove the brightest point sources and obtain 
very good amplitude estimations. The MEM technique is 
then applied to these maps to reconstruct the different com- 
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Figure 8. Input and recovered point source catalogues for the 44 GHz (top) and 545 GHz (bottom) Planck channels. The catalogues 
are convolved with the corresponding Planck beams 

ponents (except the remaining point sources contribution, proves the accuracy of the component separation of all the 

which is treated as an additional 'noise'). Following the ap- diffuse components. This is so because the MHW subtrac- 

proach discussed in H99, we generate mock observed data tion algorithm is efficient at removing the brightest point 

from our reconstructions. These maps are then subtracted sources, whereas MEM has greatly reduced the contamina- 

from the initial data. This provides data residuals maps tion due to fainter sources. We compare the reconstructions 

which mostly contain instrumental noise plus point source achieved when the MHW is or is not applied. In particular, 

emission (with slight traces of unrecovered diffuse compo- Figure |^ shows how many point sources would remain in 

nents). These residual maps are then analysed again with the reconstructed components if the MHW were not used. 

the MHW in order to refine the number of detections and We can see that a large number of point sources are removed 

the amplitude estimation of the point sources. from the dust and free-free maps. There are also a handful of 

, , , ,. , . ,. i ■ ■ point source that would contaminate the synchrotron emis- 

As already discussed m section ¥L the joint analysis un- 
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Table 4. M&M estimation of a two-thirds of the sphere point 
source catalogue (see text for details). 



Frequency 
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~ counts in 


(%) 


(GHz) 


M&Mc 


the model 




30 


5000 


7500 


66 


41 


4500 


4500 


75 


70 


5800 


6800 


85 


100 (LFI) 


7000 


10800 


65 


100 (HFI) 


12000 


24300 


49 


143 


800 


900 


89 


217 


800 


850 


94 


353 


6000 


13000 


46 


545 


20000 


32000 


63 


857 


82000 


200000 


41 



sion, coming from the low frequency channels. A lower num- 
ber of point sources would affect the CMB reconstruction 
since the cosmological signal is the main component at the 
intermediate Planck frequencies, where point source emis- 
sion is lower. Finally, a few point sources would be misiden- 
tified with SZ clusters, appearing in the reconstruction at 
the reference frequency as sharp negative features. 

In the previous section, we gave estimates of the point 
source catalogues that MEM, MHW and the joint method 
provide for these simulations (see Table ^|) . We see that the 
joint analysis provides, in general, a more complete cata- 
logue than each of the methods on its own, reaching lower 
fluxes and with point source amplitude more accurately esti- 
mated. The improvement is especially clear at high frequen- 
cies due to the high resolution of those Planck channels. The 
differences between the number of detections in the MEMc 
and M&Mc are smaller for the channels between 30 and 100 
GHz. This is due to the difficulty of detecting point sources 
using the MHW when the background has a similar scale 
variation to that of the point sources (see V01 for more de- 
tails). Hence the main contribution of the MHW at these 
frequencies is in improving the amplitude estimation. In Ta- 
ble ky we give an estimate of the number of point sources that 
would be detected with this combined method in two-thirds 
of the sky after 12 months of observation with the Planck 
satellite. This number is simply obtained by multiplying the 
counts of Table || by the ratio between the solid angle covered 
by two-thirds of the sphere and that covered by our simu- 
lations. We compare the recovered point source catalogue 
with the simulated one, with a cutoff as given by the 'Min- 
imum Flux' column for M&M given in Table |§| We can see 
that for most of them, the percentage of detection is around 
or above 50%. Current evolution models of dust emission 
in galaxies (see, e.g., Franceschini et al. 1994; Guiderdoni 
et al. 1998; Granato et al. 2000) give different predictions 
for counts in the high frequency Planck channels. On the 
other hand, all these models predict a very sharp increase of 
the far-IR/sub-mm galaxy counts at fluxes ~ 20— 100 mjy. 
Therefore, given the detection limits of Table ^, Planck data 
alone will not be able to disentangle among different mod- 
els, although it could marginally detect the sharp increase in 
the counts in the channels where the minimum flux achieved 
lies below 100 mjy. In any case, Planck will provide very 
useful data on counts, in a flux range not probed by other 



experiments. These data, complementary to the deeper sur- 
veys from the ground or from the space (ESA FIRST and 
ASTRO-F/IRIS missions), will surely allow to discriminate 
among the various evolutionary scenarios. 

Spectral information about the point sources could also 
be used to improve further the recovered catalogues. Indeed, 
V01 have shown that following point sources through adja- 
cent channels, one can estimate the spectral indices of the 
different point source populations. This would allow the re- 
covery of point sources that, albeit below the detection limit, 
have an amplitude and position in agreement with those pre- 
dicted from adjacent channels. 

Finally, as pointed out in the previous section, the flux 
limits achieved in the M&Mc are close to the noise level. 
Indeed, the faintest point sources detected in the catalogue 
have a fluxes which are 3.0, 2.4, 1.6, 1.2, 1.1, 11.2, 6.7, 1.1, 
1.3 and 1.4 times the noise rms in the 30, 44, 70, 100(LFI), 
100(HFI), 143, 217, 353, 545 and 857 GHz channels, respec- 
tively. To reach fainter fluxes in these channels is a difficult 
task, since we are very close to the noise level except for the 
143 and 217 GHz channels. On the other hand, if we sub- 
tract the MHWc sources from the original data at 143 and 
217 GHz, instead of the point sources detected by the 5a Wd 
criterion, we could greatly increase the number of sources 
and the depth of the M&Mc at those frequencies. A possi- 
bility to improve the results at all frequencies could be to 
denoise these data residual maps. One way is using wavelet 
techniques that have been proved to be very efficient at re- 
moving noise from CMB maps (Sanz et al. 1999a, b). How- 
ever, care must be taken when denoising the residual maps 
since the denoising procedure may change the profile of the 
point source in wavelet space. A detailed study of the prop- 
erties of the denoised map would then become necessary. In 
this case, instead of the Mexican Hat, one could use a cus- 
tomised pseudo-filter to detect point sources in the residual 
maps as proposed by Sanz et al (2000). 

A natural way to improve the results is to subtract the 
recovered M&Mc from the original data and applying the 
MEM algorithm again. This process could be performed 
iteratively until the flux limits and the number of counts 
converge. However, this method has some disadvantages. As 
pointed out in the previous section, if the sources subtracted 
from the input maps have a large error, this could mislead 
the MEM algorithm. This is the reason why we choose to 
subtract the catalogue achieved by the 5o~ Wd criterion in- 
stead of subtracting the one given by the 50% error criterion. 
The number of point sources with large errors in the M&Mc 
is larger than in the MHWc obtained with the 50% error cri- 
terion. Hence, a more detailed analysis becomes necessary 
in order to improve the results with an iterative approach. 
Such a study will be performed in a future work, where the 
combined technique will be extended to the sphere. More- 
over, the flux limits are already close to the noise level and 
thus we do not expect the detection levels to change sub- 
stantially (except for the 143 and 217 GHz channel, that 
can be clearly improved). 
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